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ABSTRACT 

Unified  constitutive  models  have  previously  been  developed  at  AMRL  and 
implemented  into  the  PAFEC  and  ABAQUS  Finite  Element  packages  to  predict  the 
stress-strain  response  of  structures  that  undergo  multi-axial,  rate  dependent,  cyclic 
plasticity.  In  this  work,  the  ABAQUS  version  of  the  model  has  been  enhanced  to 
account  for  advances  in  unified  constitutive  modelling,  and  for  consistency  with  the 
existing  PAFEC  version  of  the  model.  Numerous  test  cases  have  been  run  to  validate 
the  model's  implementation  in  ABAQUS.  A  full  3-dimensional  analysis  of  a  proposed 
test  specimen  has  been  undertaken  to  study  residual  stress  distribution  with  the 
ultimate  goal  of  helping  to  manage  fatigue  crack  growth  in  the  F-lll  aircraft. 
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Implementation  of  a  Unified  Constitutive 
Model  into  the  ABAQUS 
Finite  Element  Package 


Executive  Summary 

High  strength  D6ac  steel  is  used  in  the  F-lll  aircraft  for  structural  components  that 
undergo  high  stresses.  Despite  D6ac  steel's  high  strength,  some  loading  conditions  on 
the  aircraft  are  severe  enough  to  cause  localised  yielding.  Yielding  of  the  material 
leaves  residual  stresses,  which  can  lead  to  fatigue  cracking  under  normal  service 
loads.  Management  of  the  structural  integrity  of  the  F-lll  requires  knowledge  of  the 
crack  growth  rate  behaviour  at  critical  locations  in  the  aircraft.  Classical  plasticity 
modelling  techniques  have  been  shown  to  have  limitations  in  predicting  material 
response  under  non-symmetric  cyclic  loading.  Consequently,  over  the  last  decade, 
AMRL  has  researched  and  developed  a  unified  constitutive  model  to  better  predict 
the  residual  stress  field  after  a  Cold  Proof  Load  Test  (CPLT)  load  cycle. 

Unified  constitutive  models  have  previously  been  developed  at  AMRL  and 
implemented  into  the  PAFEC  and  ABAQUS  Finite  Element  packages  to  predict  the 
stress-strain  response  of  structures  that  undergo  multi-axial,  rate  dependent,  cyclic 
plasticity.  In  this  work,  the  ABAQUS  version  of  the  model  has  been  enhanced  to 
account  for  advances  in  unified  constitutive  modelling,  and  for  consistency  with  the 
existing  PAFEC  version  of  the  model.  Numerous  test  cases  have  been  run  to  validate 
the  model's  implementation  in  ABAQUS.  A  full  3-dimensional  analysis  of  a  proposed 
test  specimen  has  been  undertaken  to  study  residual  stress  distribution  with  the 
ultimate  goal  of  helping  to  manage  fatigue  crack  growth  in  the  F-lll  aircraft. 

The  result  of  the  work  reported  herein  is  the  implementation  and  validation  of  the 
best  available  material  model  for  D6ac  steel  in  the  preferred  nonlinear  Finite  Element 
package  ABAQUS.  The  successful  numerical  analysis  of  a  proposed  test  specimen 
shows  that  an  experimental  program  should  be  undertaken. 
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1.  Introduction 

Many  of  the  highly  stressed  structural  components  in  the  F-lll  aircraft  are  made  of  high 
strength  D6ac  steel.  Nonetheless,  for  some  of  these  components  some  infrequent  high 
loading  conditions  on  the  aircraft  are  sufficient  to  cause  yielding  of  the  material, 
particularly  around  stress  concentrators  such  as  holes.  The  yielding  of  the  material  leaves 
tensile  residual  stresses,  which  can  cause  fatigue  cracking  under  subsequent  normal 
service  loads.  It  has  been  found  in  both  RAAF  and  USAF1 2  service  of  the  F-lll  that 
fatigue  cracks  have  developed  at  these  yielded  locations  and  have  grown  through  the 
yield  zone  and  residual  stress  field.  The  main  loading  condition  causing  this  yielding  is 
the  Cold  Proof  Load  Test  (CPLT),  in  which  the  aircraft  is  subject  to  limit  load  cycles  in  a 
ground  test  rig.  This  test  is  performed  at  regular  intervals  and  is  the  cornerstone  of  the 
safety-by-inspection  program  for  the  D6ac  steel  components  in  the  aircraft.  CPLT 
involves  the  application  of  two  cycles  of  negative  and  positive  limit  loads  at  two  wmg 
sweep  angles  as  follows:  @  56°  wing  sweep,  0  -»  -2.4  g  -*  0  -+  +7.33  g  ->  0;  @  26°  wing 
sweep  0  — >  -3.0  g  — >  0  — »  +7.33  g  — +  0. 

Cyclic  yielding  induced  by  these  two  CPLT  load  cycles  is  complex  to  analyse.  Finite 
Element  Analysis  (FEA)  provides  a  powerful  means  of  predicting  materially  nonlinear 
response  in  complex  structures  undergoing  cyclic  loading.  However,  the  disadvantages  of 
classical  plasticity  models  used  in  FEA  to  simulate  inelastic  cyclic  response  led  AMRL  to 
investigate  other  constitutive  models.  State  variable  constitutive  models  offer  more 
accurate  predictions  of  strain  rate  and  cyclic  response  effects.  The  state-variable  based 
constitutive  model  used  at  AMRL  for  modelling  inelastic  response  of  D6ac  steel  was 
initially  developed  by  Stouffer  and  Bodner  [1]  and  Ramaswamy  [2],  and  has  subsequently 
been  jointly  enhanced  by  AMRL,  The  University  of  Melbourne  and  Monash  University3. 
The  model,  referred  to  as  the  unified  constitutive  model,  has  been  implemented  into  the 
PAFEC  Finite  Element  package  [3],  which  was  the  only  relevant  package  available  at  the 

time. 

Whilst  PAFEC  has  facilities  for  nonlinear  analysis,  and  will  allow  the  installation  of  user 
material  models,  the  ABAQUS  Finite  Element  package  has  developed  a  wide  reputation 
as  the  preferred  nonlinear  finite  element  code.  Indeed,  MacNeal-Schwendler  Corporation, 
the  developers  of  NASTRAN,  markets  and  recommends  ABAQUS  for  nonlinear  analysis. 
Also,  the  current  PAFEC  implementation  of  the  unified  constitutive  model  requires 
modifying  more  than  20  existing  PAFEC  routines  as  well  as  creating  more  than  30  new 
routines.  Hence,  limitations  of  the  PAFEC  implementation  are  the  time  and  expertise 
required  maintaining  the  unified  constitutive  model  when  new  versions  of  PAFEC  are 
released.  Therefore,  when  ABAQUS  became  available,  an  implementation  of  the  unified 
constitutive  model  into  that  package  was  desirable. 


1  Royal  Australian  Air  Force. 

2  United  States  Air  Force. 

3  Reference  [3]  contains  more  information  on  the  history  of  the  model's  development. 
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This  report  describes  the  implementation  of  the  Ramaswamy-Stouffer-Bodner  unified 
constitutive  model  into  ABAQUS,  which  was  based  on  development  work  done  by 
Monash  University  [4].  Section  2  discusses  limitations  of  existing  ABAQUS  plasticity 
models,  gives  the  mathematical  background  of  the  unified  constitutive  model,  and 
describes  the  implementation  of  the  model  into  ABAQUS,  noting  differences  from 
previous  work  [4],  Section  3  describes  a  proposed  test  specimen,  and  details  preliminary 
(elastic)  finite  element  modelling  of  the  specimen.  In  Section  4  an  elasto-plastic  analysis  of 
the  specimen  is  reported,  with  the  aim  of  characterising  the  fatigue  crack  growth  through 
a  residual  stress  field  left  after  yield  of  a  D6ac  component. 


2.  Constitutive  Modelling  Approach 


2.1  Existing  ABAQUS  Plasticity  Models 

ABAQUS  has  several  options  for  plasticity  models  [5]  that  are  applicable  to  the  inelastic 
behaviour  of  metals.  Most  of  these  models  are  incremental,  in  that  the  total  strain  is 
separated  into  an  elastic  part,  on  which  all  changes  in  stress  depend,  and  a  plastic  part. 
Incremental  plasticity  models  generally  are  defined  by  the  following:  (i)  a  yield  function, 
(ii)  a  flow  rule,  and  (iii)  a  hardening  law. 

ABAQUS  has  two  yield  criteria  suitable  for  predicting  behaviour  of  metals,  viz.,  (i)  Von 
Mises  criterion,  and  (ii)  Hill  criterion,  which  is  a  form  of  the  Von  Mises  function  that 
allows  for  anisotropic  response.  The  ABAQUS  material  models  for  plasticity  adopt 
associative  flow  rules,  viz.  plastic  straining  is  normal  to  the  yield  surface.  ABAQUS 
allows  for  perfectly  plastic  response,  isotropic  hardening  and  kinematic  hardening,  and 
has  a  nonlinear  isotropic/kinematic  hardening  model  which  is  intended  to  model  the 
cyclic  loading  of  metals.  Strain  rate  effects  can  be  accounted  for  in  ABAQUS' s  simple 
isotropic  hardening  scheme,  but  are  not  available  for  kinematic  hardening  and  combined 
isotropic/ kinematic  hardening.4  Creep  models  are  available  in  ABAQUS,  but  these  are 
not  applicable  to  creep  under  cyclic  loading.  Also,  if  creep  and  plasticity  occur 
simultaneously,  the  response  may  not  be  predicted  correctly. 

Plasticity  in  ABAQUS  is  one  of  a  range  of  nonlinear  material  behaviours  that  can  be 
simulated.  The  general  approach  in  ABAQUS  to  the  solution  of  nonlinear  problems  is  to 
apply  the  loading  in  steps,  with  the  load  in  each  step  being  divided  into  increments. 
Equations  representing  the  response  of  the  structure  to  a  load  increment  are  solved  by 
iteration  -  the  Newton-Raphson  method  (and  variations)  is  adopted,  with  judicious 
monitoring  of  progress  towards  convergence  to  ensure  a  computationally  efficient 


4  When  strain  rate  effects  are  negligible,  ABAQUS's  combined  isotropic/ kinematic  hardening 
(nonlinear  kinematic  hardening)  model  may  provide  an  efficient  alternative  to  the  Unified 
Constitutive  Model. 
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solution.  After  a  material  point  has  yielded,  its  stress  components  are  updated  for  the 
increment  using  an  implicit  integration  scheme. 

A  powerful  feature  of  the  ABAQUS  Finite  Element  Package  is  the  ability  to  include  a 
users  constitutive  model  through  the  UMAT  facility. 


2.2  The  Unified  Constitutive  Model  -  Background 


As  shown  in  the  previous  section,  ABAQUS  plasticity  modelling  has  several 
disadvantages,  in  common  with  all  classical  techniques.  Unified  constitutive  models  offer 
improved  representation  of  strain  rate  and  cyclic  response  effects  [6],  which  means  time 
dependent  responses  such  as  creep,  stress  relaxation  and  cyclic  hardening/ softening  can 
be  predicted  with  greater  accuracy. 

Unified  constitutive  models5  are  so-called  because  classical  plastic  strain  and  creep  strain 
are  'unified'  into  a  single  inelastic  strain  s\. .  This,  by  definition,  overcomes  the  problems 

of  classical  plasticity  when  plasticity  and  creep  occur  together.  Also,  with  no  distinction 
between  plasticity  and  creep  regimes,  high  and  low  temperature  responses  are  simulated 
with  the  same  set  of  basic  equations.  Thus,  in  isothermal  applications,  the  total  strain  ( £- ) 


is, 


£,  -  <+< 


(1) 


where  sf.  is  the  elastic  strain.  The  total  stress  ( cr. )  is  defined  by, 

O'.  =  C 

ij  tjkl  ij 

where  CijkJis  the  (constant)  elastic  constitutive  tensor.  In  the  displacement-based  finite 

element  method,  total  strains  are  available  at  each  stage  of  the  solution.  The  major 
function  of  a  unified  constitutive  model  is  to  predict  the  values  of  inelastic  strain,  from 
which  elastic  strain  and  total  stress  can  be  determined.  Inelastic  strain  is  modelled  (in  rate 

form)  as  a  function  of  deviatoric  stress  (Sy)  and  two  state  variables,  back  stress  (Cl  ), 


and  drag  stress  (Z) .  Back  stress  is  a  tensor  quantity  that  represents,  at  microscopic  level, 
dislocation  movements  and  pile-ups.  Drag  stress  is  a  scalar  that  models  the  evolution  of 
subgrains  and  other  orientation-independent  microscopic  properties.  The  inelastic  flow 
equation  of  the  unified  constitutive  model  developed  initially  by  Stouffer  and  Bodner  [1] 
and  Ramaswamy  [2]  is: 

s‘.  =  Dexp 


1 

r  z>  Y 

to-n.) 

2 

Uk2J 

VkT 

5  The  background  to  the  unified  constitutive  model  detailed  in  this  section  follows  the  description 
given  by  Trippit  [7]. 
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where  £-  is  the  inelastic  strain  rate,  D  and  n  are  temperature  and  strain  rate  dependent 
material  parameters,  and  K2is  the  second  invariant  of  the  over  stress  tensor 
(Os  =  S0.  -Q..),viz., 

k2  =  ^s,-q,us,-Q,)  « 

The  deviatoric  stress  (Sy )  is  defined  by, 

S  -  a..-  8  a 

ij  ij  V  m 

where  8;j  is  the  Kronecker  delta,  and  (Jm  is  the  mean  (or  hydrostatic)  stress,  viz.. 


1 


CTm  =  ~(<7U  +^22  +  (73}) 
3 


(6) 


In  Equation  (3),  for  D6ac  steel  at  room  temperature,  D  =  10  000  for  tensile  and  fatigue 
loadings  (strain  rates  up  to  1.0  sec  -*),  and  n  =  3  to  model  the  strain  rate  insensitivity  of  the 
response.  The  inelastic  flow  equation  (Eqn  3)  requires  values  of  back  stress.  Several  forms 
of  the  back  stress  evolution  equation  have  been  proposed  for  the  unified  constitutive 

model.  In  this  development  the  back  stress  rate  (Cly)  is  defined  by  the  multi-axial  stress 
form  of  an  equation  derived  by  Trippit  [7],  p.79, 

£V  a 


where  ft  is  a  hardness  related  material  parameter,  Q„is  the  steady  state  value  of  back 
stress6,  f2  is  another  material  parameter,  is  the  equivalent  inelastic  strain  rate  and 
S.,  is  the  deviatoric  stress  rate.  The  equivalent  inelastic  strain  rate  is  defined  as  : 


£equ  = 


2 

'3  11  1 


(8) 


V 

Also  required  for  the  evaluation  of  inelastic  strain  rates  is  the  drag  stress  (Z).  Its  evolution 
equation  is : 

Z  =  Z  j  +  (Z  0  -  Z  j  )exp(-m£,e/qu )  (9) 

where  Z0  and  Z,  are  initial  and  final  values  of  drag  stress,  m  controls  the  rate  at  which 
Z  evolves  from  Z0toZ,,  and  £  qu  is  the  equivalent  inelastic  strain.  Z0,  Z,  and  m  are 

material  parameters,  and,  from  Equation  (8),  £equ  is  given  by: 


<u = 


(10) 


6  Q  is  also  referred  to  as  'equilibrium  back  stress'  or  'saturated  back  stress  [6].  For  tensile  and 
fatigue  response,  =  Qmax ,  the  maximum  back  stress  reached  at  inelastic  saturation. 
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Equations  (3),  (7)  and  (8)  form  a  system  of  ordinary  differential  equations  (ODEs),  the 
solutions  of  which  yield  components  of  inelastic  strain  and  back  stress,  as  well  as  the 
equivalent  inelastic  strain.  The  Livermore  Solver  [8]  has  proved  to  be  very  effective  in  the 
solution  of  the  unified  constitutive  model  ODEs  [7],  and  is  used  in  this  implementation  of 
the  model.  Searl  and  Paul  [9]  have  derived  values  for  the  material  parameters 

f  f  ,  Q  (=  Q  ),  Zn ,  Z.  and  m  for  D6ac  steel. 

J'l  ’  2 »  max'.  «/’  0>  1 


2.3  The  Unified  Constitutive  Model  -  Implementation 

The  unified  constitutive  model  has  previously  been  implemented  at  AMRL  in  the  PAFEC 
finite  element  package  [3].  That  implementation  is  part  of  work  on  the  Ramaswamy- 
Bodner-Stouffer  version  of  the  unified  constitutive  model  that  has  been  carried  out  at 
AMRL  and  elsewhere  over  the  past  ten  years7.  The  PAFEC  implementation  of  the  model 
for  D6ac  steel  involved  new  routines  for  automatic  time  step  load  calculation  and 
structural  convergence  control,  as  well  as  routines  to  model  material  behaviour. 
However,  ABAQUS  is  able  to  completely  manage  the  solution  process  at  the  structural 
level.  Hence,  this  implementation  only  required  the  development  of  UMAT  (User 
Material)  modules  to  calculate  total  stress  and  a  consistent  Jacobean  matrix  for  a  material 
point.  Monash  University  was  contracted  to  implement  the  unified  constitutive  model 
developed  at  AMRL  into  ABAQUS  [4].  The  implementation  presented  herein  is  based  on 
that  work. 


Several  modifications  have  been  made  to.  the  Monash  version  for  the  current  ABAQUS 
installation,  the  main  purposes  being  to  update  the  model  in  accord  with  advances  in 
unified  constitutive  modelling,  and  to  produce  a  material  model  which  will  give  results 
consistent  with  equivalent  PAFEC  analyses.  The  back  stress  evolution  equation  has  been 
updated.  The  previously  used  one  was  (compare  to  Eqn  7): 


Q .  =  f,£- ---Ufi.fi,,  +  f,5, 


2  0. 


equ 


2  y 


(ii) 


A  measure  of  effective  inelastic  strain  is  required  in  the  drag  stress  evolution  equation 
(Eqn  9).  The  previous  implementation  used  a  different  definition  of  effective  inelastic 

strain,  viz.,  the  inelastic  strain  magnitude  ( ),  which  is  defined  by, 


(12) 


The  drag  stress  should  evolve  during  loading  by  asymptoting  to  a  saturated  value. 
However,  during  cyclic  loading  the  inelastic  strain  magnitude  will  increase,  then  decrease 
and  increase  again,  resulting  in  fluctuating  values  of  drag  stress.  The  previous  formulation 
prevented  the  fluctuation  by  modifying  the  application  of  the  drag  stress  equation. 
However,  equivalent  plastic  strain  (Eqn  10),  which  is  used  in  the  current  PAFEC  (L8.1) 
code,  monotonically  increases  during  loading.  Adopting  this  definition  of  equivalent 
inelastic  strain  results  in  a  more  general,  more  efficient  and  simpler-to-use  formulation. 


7  Searl  and  Paul's  report  [9]  contains  a  comprehensive  list  of  references  on  the  development  of  the 
model  and  its  implementation  into  PAFEC  at  AMRL. 
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because  strain  softening  can  be  modelled,  fewer  state  variables  are  needed,  and  Initial 
Conditions  are  no  longer  needed  in  the  ABAQUS  input  deck  (for  tension  and  fatigue 
loadings)  Other  modifications  are  an  increase  in  the  number  of  ABAQUS  element  types 
that  can  use  the  material  model  (Section  2.4),  and  the  facility  to  analyse  a  model 
containing  more  than  one  element  type. 

2.4  Supported  Element  Types 

The  ABAQUS  user  material  interface  provides  the  numbers  of  direct  stresses  and  shear 
stresses  as  input  to  the  user-written  constitutive  model.  In  the  most  general  case,  there  are 
three  direct  stresses  and  three  independent  shear  stresses  at  a  material  point.  However, 
most  element  types  do  not  require  the  full  set.  Hence  the  constitutive  model  should  be 
developed  for  variations  in  the  number  of  stresses.  This  implementation  of  the  unified 
constitutive  model  supports  all  Version  5.6  ABAQUS  elements  with  displacement 
degrees-of-freedom,  as  summarised  in  Table  1. 

2.5  Validation 

2.5.1  Introduction 

It  is  necessary  to  verify  that  the  implementation  of  the  unified  constitutive  model  into 
ABAQUS  is  performing  as  expected.  One,  two  and  three-dimensional  test  cases  have  been 
created  to  prove  the  implementation  -  these  tests  are  similar  to  the  suite  developed  to 
validate  the  installation  of  the  model  into  PAFEC  [3].  In  both  instances  the  test  case 
results  have  been  compared  to  output  from  a  standalone  FORTRAN  program  which 
simulates  uniaxial  response.  A  description  of  this  program,  which  is  named  TREACLE,  is 
given  in  Appendix  A.  The  application  of  the  unified  constitutive  model  to  D6ac  steel  has 
been  validated  by  Searl  and  Paul  [9],  where  appropriate  material  constants  have  been 
determined,  and  experimental  test  data  and  numerical  results  compared. 

2.5.2  Description  of  Test  Cases 

The  test  cases  used  element  types  that  encompassed  all  relevant  variations  of  direct  and 
shear  stresses,  as  outlined  in  Table  1.  Tests  also  involved  a  variety  of  loadings  and 
different  configurations  of  elements.  The  loading  types  were  uniaxial  pull,  constrained 
pull,  shear  and  uniaxial  cyclic.  Element  shapes  were  cubic,  square  or  linear,  with  a 
constant  element  edge  length  of  1  mm. 
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Table  1.  Summary  of  Supported  Element  Types 


Element  type 

Number  of  direct  stresses 

Number  of  independent 
shear  stresses 

3D  Solid 

3 

3 

Plane  Strain 

3 

1 

Plane  Stress,  Shells, 
Axisymmetric  shells  with 
non-axisymmetric  loading. 
Elbows 

2 

1 

Axisymmetric  shells 

2 

0 

Beams  with  torsional  shear 
stress 

1 

1 

Truss,  Beams  without 
torsional  shear  stress 

1 

0 

Tables  2,  3  and  4  give  a  summary  of  these  test  cases  with  the  types  of  elements  and 
loading  conditions  that  have  been  validated.  In  test  cases  where  the  response  was 
uniaxial  and  uniform  the  solutions  were  compared  to  answers  from  the  TREACLE 
program  In  all  tests,  criteria  involving  stress  and  strain  values  were  invoked  for 
validating  the  test  results.  In  Tables  2  and  3  the  test  case  numbers  1-13  perform  the  same 
verifications  as  the  corresponding  PAFEC  tests  [3].  Cases  14-17  are  further  tests.  Material 
constants  adopted  in  the  test  cases  were  the  same  as  those  used  in  the  PAFEC  validation 

[3],  viz., 

Qmax  =1250  MPa  Z0  =250  MPa 

m  =9  Z\  =  550  MPa 

E  =  207  GPa  fi  =72 


=  0.32 


f2  =  0.68 


The  maximum  inelastic  strain  rate  (D)  and  strain  rate  sensitivity  parameter  (n)  have  the 
values  defined  in  Section  2.2.  Figure  1  is  the  typical  stress-strain  curve  for  the 
aforementioned  data,  to  40  000  ps  at  a  strain  rate  of  100  gs  sec  L  The  data  points  of  which 
were  generated  using  TREACLE ,  which  is  described  in  Appendix  A. 
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Table  2.  Summary  of  Test  Cases  with  3-Dimensional  Elements 


Test 

no. 

Data 

file 

Element 

type 

Element 

configu¬ 

ration 

Loading 

condition 

Description 

1 

A3D1 

C3D20R 

2x2x2 

Free  Uniaxial 

Test  the  general  stress  to  strain 

Pull 

response. 

2 

A3D2 

C3D20R 

2x2x2 

Constrained 

Test  the  general  stress  to  strain 

Pull 

response  under  different 

hydrostatic  conditions  than  test  1. 

3 

A3D3 

C3D8R 

1 

Shear 

Test  the  general  stress  to  strain 

response  in  shear. 

4 

A3D4 

C3D20R 

2x2x2 

Free  Uniaxial 

Test  the  cyclic  stress  to  strain 

Cyclic 

response. 

Table  3. 

Summanj  of  Test  Cases  with  2-Dimensional  Elements 

Test 

Data 

Element 

Element 

Loading 

Description 

no. 

file 

type 

configu- 

condition 

ration 

5 

APS1 

CPS8 

2x2 

Free  Uniaxial 

Plane  stress.  Test  the  general  stress 

Pull 

to  strain  response. 

6 

APS2 

CPS8 

2x2 

Constrained 

Plane  stress.  Test  the  general  stress 

Pull 

to  strain  response  under  different 

hydrostatic  conditions  than  test  5. 

7 

APS3 

CPS4 

1 

Shear 

Plane  stress.  Test  the  general  stress 

to  strain  response  in  shear. 

8 

APS4 

CPS8 

2x2 

Free  Uniaxial 

Plane  stress.  Test  the  cyclic  stress 

Cyclic 

to  strain  response. 

9 

APE1 

CPE8 

2x2 

Free  Uniaxial 

Plane  strain.  Test  the  general  stress 

PuH 

to  strain  response. 

10 

APE2 

CPE4 

1 

Shear 

Plane  strain.  Test  the  general  stress 

to  strain  response  in  shear. 

11 

APE3 

CPE8 

2x2 

Free  Uniaxial 

Plane  strain.  Test  the  cyclic  stress 

Cyclic 

to  strain  response. 

12 

ASH1 

S8R 

2x2 

Free  Uniaxial 

Test  the  general  stress  to  strain 

PuH 

response. 

13 

AMU1 

S8R 

l+l+l8 

Free  Uniaxial 

Multiple  element  type  tests. 

CPS8 

Pun 

CPE8 

14 

ASA1 

SAX1 

1 

Free  Uniaxial 

Test  the  general  stress  to  strain 

Pun 

response. 

8  One  element  of  each  type,  not  connected  to  each  other,  individually  loaded. 
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Table  4.  Summary  of  Test  Cases  with  1-Dimensional  Elements 


Test 

no. 

Data 

file 

Element 

type 

Element 

configu¬ 

ration 

Loading 

condition 

Description 

15 

AB1 

B31 

1 

Free  Uniaxial 
Pull 

Test  the  general  stress  to  strain 
response. 

16 

AB2 

B31 

1 

Torsion 

Test  the  general  stress  to  strain 
response  in  pure  torsional 
loading. 

17 

ATI 

T2D2 

1 

Free  Uniaxial 
Pull 

Test  the  general  stress  to  strain 
response. 

Two  of  the  test  cases  are  described  below.  The  model  for  Test  1,  which  is  summarised  in 
Table  2,  was  a  three  dimensional  block  composed  of  eight  3D  continuum  elements.  The 
block  was  loaded  by  prescribed  displacements  on  four  of  the  element  faces,  as  shown  in 
Figure  2a.  Figure  2b  shows  the  block  restrained  on  the  faces  opposite  to  the  loaded  faces, 
as  well  as  on  two  adjoining  planes,  to  prevent  rigid  body  motion  and  preserve  uniformity 
of  response.  The  applied  loading  caused  a  maximum  strain  of  40  000  ps.  The  time  of  the 
analysis  was  400  sec,  giving  a  constant  strain  rate  of  100  ps  sec  -k  Results  of  this  analysis 
are  shown  in  Figure  3.  Test  4,  which  is  summarised  in  Table  1,  had  the  same  finite 
element  configuration,  maximum  prescribed  displacement  and  restraints  as  the  Test  1, 
but  underwent  the  cyclic  loading  shown  in  Figure  4.  Results  of  this  analysis  are  shown  in 
Figure  5.  Results  of  all  tests  satisfied  criteria  for  stress  and  strain,  and,  where  applicable, 
results  agreed  with  output  from  the  standalone  TREACLE  program,  as  shown,  for 
example,  in  Figures  3  and  5. 

3.  Elastic  Numerical  Analysis  of  Test  Specimen 


Elastic  and  inelastic  numerical  analyses  in  this  study  have  been  performed  using  the 
ABAQUS  Finite  Element  Package,  Version  5.6,  running  on  a  Hewlett-Packard  K  Series 
9000  computer  at  AMRL.  Patran  Version  7.5  has  been  used  for  all  pre-  and  post¬ 
processing.  Twenty-noded  hexahedral  elements  have  been  used  for  all  runs. 

3.1  Geometry  of  Test  Specimen 

Figure  6  shows  the  geometry  of  the  test  specimen.  The  specimen  is  made  from  D6ac  steel 
bar  stock  and  has  outside  dimensions  550  by  150  by  30  mm.  The  middle  section  has  been 
reduced  in  thickness  to  5  mm,  and  a  30  mm  diameter  hole  has  been  drilled  in  the  centre  of 
the  plate  to  act  as  a  stress  concentrator.  Section  A-A  in  Figure  6  shows  that  30  mm  wide 
'flanges'  remain  at  the  edge  of  the  central  plate,  to  give  the  specimen  stability  when  tested 
in  compression. 
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3.2  Preliminary  Analyses 

3.2.1  Introduction 

Several  preliminary  elastic  analyses  were  conducted.  The  aims  of  these  analyses  were  to 
determine  the  following:  (i)  the  most  appropriate  way  to  model  the  load,  (ii)  if  lateral 
stability  would  be  a  limiting  factor  when  testing  in  compression,  and  (ui)  what  was  likely 
to  be  an  adequate  mesh  density  for  detailed  analyses.  Because  none  of  these  aims 
involved  an  accurate  prediction  of  the  stresses  around  the  hole,  the  analyses  could  be 
carried  out  with  a  relatively  coarse  mesh.  However,  a  refined  mesh  model  would  also  be 
needed  to  obtain  accurate  estimates  of  stress  for  elastic  and  inelastic  cases.  To  facilitate 
general  refinement  of  the  mesh  and  to  give  flexibility  in  local  mesh  refinement  a  blocking 
scheme  was  used.  Figure  7  shows  the  blocking  scheme  for  a  1/8-model  of  the  test 

specimen. 

3.2.2  Modelling  of  Load  Application 

For  investigating  the  load  application,  symmetry  allows  a  1/8-model  to  be  used.  Figure  8 
shows  the  model,  which  has  495  20-noded  hexahedral  elements,  with  Young  s  Modulus 
equal  to  210  GPa,  and  Poisson's  Ratio  0.3.  The  restraints  introduced  by  modelling  /  o 
the  test  specimen  are  sufficient  to  constrain  the  model  against  rigid  body  motion. 

The  aims  of  this  set  of  analyses  were  to  determine  the  following:  (i)  whether  the  size  of 
the  finite  element  model  could  be  reduced  by  neglecting  the  part  of  the  specimen  that 
would  be  gripped  in  the  testing  machine,  (ii)  if  the  contact  pressure  applied  through  the 
grips  of  the  testing  machines  affected  stresses  at  the  hole,  and  (iii)  if  the  assumed 
distribution  of  the  tensile/ compressive  force  affected  stresses  at  the  hole. 

The  following  five  load  cases  were  analysed: 

i.  Base  case  -  load  applied  to  the  end  surface  (150  mm  by  30  mm  surface,  refer  to 

Figure  6)  as  a  pressure. 

ii.  Short  sample  -  the  150  mm  by  140  mm  by  30  mm  (refer  to  Figure  6)  portions  at  the 
ends  of  the  sample  are  removed.  The  load  is  applied  as  an  end  pressure  ,  as  in  the 

base  case. 

iii.  Contact  pressure  -  load  applied  to  the  end  surface  as  a  pressure,  with  the  testing 
machine  grips  exerting  a  uniform  contact  pressure  of  30  MPa. 

iv.  Uniform  Shear  -  load  applied  to  the  surfaces  of  the  specimen  in  contact  with  the 
testing  machine  grips,  as  a  uniformly  distributed  shear  force. 

v.  Linear  Shear  -  load  applied  to  the  surfaces  of  the  specimen  in  contact  with  the 
testing  machine  grips,  as  a  shear  force  with  a  bilinear  distribution,  viz.,  at  a 
maximum  half  way  along  each  grip. 
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In  all  cases,  the  equivalent  longitudinal  load  was  135  kN  (or  30  MPa).  Table  5  shows 
values  of  maximum  stress  for  each  load  case9,  and  the  difference  from  the  base  result. 


Table  5  -  Results  of  Load  Application  Modelling 


Load  case 

Maximum  stress  (MPa) 

%  difference  (from  (1)) 

1 

369.3 

- 

2 

415.4 

+12.5 

3 

372.0 

+0.7 

4 

369.6 

+0.1 

5 

369.6 

+0.1 

Specific  conclusions  from  these  analyses  are  the  following: 

i.  the  'short  sample'  results  differ  significantly  from  the  base  case.  Hence,  a  full-length 

model  must  be  used. 

ii.  addition  of  contact  pressure  has  a  minor  effect  on  results  (about  0.7%).  Therefore,  it 
can  be  disregarded  in  modelling. 

iii.  change  in  application  of  load  from  end  pressure  to  shear  force  makes  little  difference 

(less  than  0.10%).  Changing  the  shear  force  distribution  from  uniform  to  bilinear 
makes  no  discernible  difference. 

The  general  conclusion  is  that  the  test  specimen  should  be  modelled  to  full  length,  with 
the  load  applied  as  an  end  pressure,  and  no  contact  pressure  exerted  by  the  grips. 


3.2.3  Buckling 

The  first  buckling  analysis  was  done  with  a  full  model  (3960  elements),  so  that  no 
buckling  modes  would  be  excluded  by  assumptions  of  symmetry.  Young's  Modulus  was 
chosen  as  209  GPa,  and  Poisson's  Ratio  as  0.3.  The  critical  load  that  the  test  specimen 
could  reach  was  considered  to  depend  greatly  on  the  degree  of  restraint  provided  by  the 
jaws  of  the  test  machine.  It  was  decided  to  initially  adopt  a  simple  (conservative) 
assumption  for  the  degree  of  restraint,  and  improve  modelling  of  the  restraint  only  if  the 
results  were  unsatisfactory.  The  first  condition  was  a  hinge  at  both  ends  of  the  specimen. 
The  buckling  load  predicted  by  ABAQUS  in  this  case  was  376  kN.  Because  the  proposed 
loads  were  to  exceed  1000  kN,  this  was  not  satisfactory. 

The  lowest  three  buckling  modes  were  examined  for  the  'hinged-at-end'  case,  and  they 
showed  that  further  analysis  could  reliably  be  undertaken  with  a  Vi-model  of  the 
specimen.  This  model  (which  was  also  modified  in  a  minor  way)  contained  1032 
elements.  After  confirming  that  the  Vi-model  gave  the  same  results  as  the  full  model  for 
the  'hinged-at-end'  case,  a  second  analysis  was  done  with  restraints  modelled  as  hinges 
halfway  along  each  grip.  The  buckling  load  in  this  case  was  589  kN,  which  was  still  not 


9  The  maximum  stresses  are  predicted  at  the  edge  of  the  hole,  nearest  the  specimen  'flanges ,  at 
half-thickness  of  the  plate.  They  act  in  the  longitudinal  direction  of  the  specimen. 
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enough.  Hence,  it  was  necessary  to  model  the  degree  of  restraint  more  accurate  y. 
Discussions  with  the  test  engineer  concluded  that  all  of  the  ends  of  the  specimen  would 
be  held  in  the  jaws  of  the  test  machine,  and  that  its  grips  would  fix  the  specimen  agarns 
all  rotations  as  well  as  translation.  Hie  buckling  load  predicted  for  this  hilly  fvced l  option 
was  2065  kN.  Figure  9  shows  that  the  deformation  mode  is  buckling  of  the  5  mm  thick 
central  section.  This  analysis  shows  that  the  test  specimen  (given  correct  alignment)  could 
be  loaded  to  approximately  1000  kN  without  the  likelihood  of  excessive  lateral 

displacements. 

3.2.4  Mesh  Density 

Figure  10  shows  the  'maximum  differences'  in  the  longitudinal  stresses  for  the  base  case, 
in§the  region  of  the  hole.  'Maximum  differences'  are  calculated  by  extrapolating 
longitudinal  stresses  in  an  element  to  its  nodes,  and  comparing  the  values  at  each  node 
from  ah  elements  which  'use'  that  node.  Hence,  it  is  a  measure  of  the  coarseness  of  the 
mesh,  and  can  be  interpreted  as  an  indication  of  the  mesh  error. 

The  size  of  the  mesh  errors  showed  that  the  mesh  needed  some  refining.  Five  elements 
were  used  through  the  half-thickness  in  the  coarse  mesh.  Because  of  the  relatively  smaU 
variation  in  longitudinal  stress  through  the  thickness,  five  elements  were  considered  to  be 
adequate  In  the  coarse  mesh  the  element  densities  in  the  two  blocks  surrounding  the  ho  e 
Con  the  1/8-model)  were  3  by  3  and  2  by  3.  For  the  refined  mesh  it  was  decided  to  try  9  by 
8  and  6  by  8.  In  general,  element  side  lengths  vary  within  a  block.  The  element  sizing  was 
based  on  a  qualitative  assessment  of  stress  gradients  for  Load  Case  1,  with,  of  course 
shorter  edge  lengths  where  stress  gradients  are  high.  The  refmed  mesh  1/8-  model  of 
test  specimen  contains  1854  elements. 

3.3  Analyses  with  Refined  Mesh 

The  base  case  in  Section  3.1.3  was  analysed  again  with  the  refined  mesh  shown  m  Figure 
11  The  maximum  stress  reached  was  363.9  MPa  (compared  to  369.3  MPa  m  die  coarse 
mesh  analysis,  a  1.5%  difference),  although  the  variation  was  greater  m  other  areas. 
Figure  12  shows  the  'maximum  differences'  in  the  longitudinal  stresses  m  the  region 
around  the  hole.  The  highest  values  in  the  area  of  interest  are  about  2  MPa,  from  whic  l 
can  be  concluded  that  the  mesh  density  is  adequate.  To  validate  conclusions  that  were 
based  on  the  coarse  mesh  (Section  3.1.3),  an  analysis  was  done  with  a  combination  of  a 
uniform  shear  force  and  a  contact  pressure  (6  MPa)  usmg  the  refmed  mesh.  The  resul 
was  a  maximum  stress  of  365.5  MPa,  only  0.4%  greater  than  the  base  case,  which  confirms 
the  conclusions  reached  in  Section  3.1.3. 
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4.  Elasto-Plastic  Numerical  Analysis  of  Test  Specimen 

4.1  Analysis  Using  Existing  ABAQUS  Plasticity 

One  of  the  main  steps  in  the  proposed  experimental  program  is  to  load  the  specimen  so 
that  a  maximum  compressive  strain  of  20  000  pe  is  reached.  Hence,  this  strain  level  must 
be  attained  in  the  numerical  analysis  with  the  unified  constitutive  model.  The  classical 
plasticity  models  available  in  ABAQUS  are  algorithmically  simple  in  comparison  with  the 
unified  constitutive  model,  and  therefore  analysis  using  a  classical  plasticity  model  is 
computationally  more  efficient.  Hence,  it  was  decided  that  an  analysis  with  classical 
plasticity  would  be  valuable  in  estimating  the  load  (and  its  incrementation)  to  be  applied 
in  a  unified  constitutive  model  analysis.  Also,  the  results  would  provide  a  limited 
validation  of  responses  predicted  by  the  unified  constitutive  model,  particularly  up  to 
initial  yield. 

The  fine  finite  element  model  used  for  elastic  analysis  was  adopted,  with  reduced 
integration  elements10.  The  plasticity  model  applied  is  the  Von  Mises  yield  criterion,  with 
Prandl-Reuss  flow  equations,  and  no  hardening  (viz.,  perfectly  plastic  behaviour). 
Young's  Modulus  (E)  is  203  GPa,  Poisson's  Ratio  (v)  is  0.3  and  the  yield  stress  is  1460 
MPa,  which  equates  to  a  yield  strain  of  approximately  7200  ps.  Figure  13  shows  the  total 
longitudinal  strains  at  a  load  of  1012  kN  (in  compression)  predicted  by  the  classical 
analysis.  The  total  load  (1350  kN)  was  applied  in  two  steps,  the  first  being  one  increment 
of  60%  of  the  total  load,  the  second  of  eight  increments  of  the  remaining  40%.  Because 
D6ac  steel  shows  significant  strain  hardening  and  this  analysis  assumes  perfectly  plastic 
response,  the  strains  predicted  at  the  load  level  shown  are  expected  to  be  higher  than 
those  attained  in  a  test.  However,  it  is  estimated  that  the  desired  strain  (20  000  ps)  would 
be  attained  at  a  load  of  approximately  1000  kN. 

4.2  Analyses  Using  Unified  Constitutive  Model 

4.2.1  Introduction 

The  finite  element  mesh  and  element  type  for  analyses  using  the  unified  constitutive 
model  were  the  same  as  for  classical  plasticity  analyses  (Section  4.1).  The  material 
constants  and  the  loading  sequence  are  described  in  the  following  sections. 

4.2.2  Constitutive  Constants 

The  material  constants  adopted  in  these  analyses  are  based  on  the  constitutive  constants 
derived  by  Searl  and  Paul  [9,  Section  7.2]  and  are  as  follows: 

Qmax  =1250  MPa  Z0  =1081  MPa 


10  The  fully  integrated  form  of  the  ABAQUS  20-noded  hexahedral  element  exhibits  volumetric 
locking  when  response  is  incompressible,  as  occurs  for  the  plasticity  model  adopted  here. 
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m  =  120  Zi  =  600  MPa 

E  =  205  GPa  fi  = 135 

v  =  0.32  f2  =  0.585 

Note  that  the  implementation  of  the  unified  constitutive  model  into  ABAQUS  allows  fi 
and  f2  to  evolve  in  a  similar  manner  to  the  drag  stress  (Eqn  4).  However,  here  the 
recommended  values  for  fi  and  f2  are  constant.  This  can  be  accomplished  by  setting  the 
initial  and  final  values  of  fi  and  f2  to  their  constant  values,  in  which  case  die  parameters 
that  would  control  the  evolution  of  fi  and  f2are  arbitrary.  As  mentioned  in  Section  2.2  the 
strain  rate  sensitivity  parameter  (n)  is  set  to  3.0.  Finally,  the  Jacobean  factor,  which 
controls  the  accuracy  of  the  Jacobean  matrix  calculated  by  the  model,  is  set  to  3. 


4.2.3  Loading  Scheme 

The  loading  scheme  for  numerical  analysis  was  intended  to  simulate  the  loading  in  the 
Plastic  Strain  Survey  of  the  proposed  CPLT  on  the  specimen.  To  this  end,  the  model  was 
first  analysed  in  monotonic  compression  to  determine  the  load  which  would  cause  a 
maximum  direct  strain  of  20  000  ps.  The  result  was  a  load  of  1012  kN.  The  loading  scheme 
for  the  numerical  analysis  is  as  follows,  with  ten  equal  increments  in  all  load  steps,  except 
die  fourth,  which  has  a  single  increment,  because  of  the  relatively  small  increase  from  9 

to  1012  kN. 


Sequence  of  loading  is  as  follows: 

i.  Ramp  the  load  in  tension  to  310  kN. 

ii.  Unload  to  zero  load. 

iii.  Ramp  the  load  in  compression  to  950  kN. 

iv.  Ramp  the  load  in  compression  to  1012  kN. 

v.  Unload  to  zero  load. 

vi.  Ramp  the  load  in  tension  to  414  kN. 

vii.  Unload  to  zero  load. 

viii.  Ramp  the  load  in  compression  to  1012  kN. 

ix.  Unload  to  zero  load. 

The  loading  rate  is  a  uniform  4.5  kN/ sec. 


4.2.4  Results 

Figure  14  shows  total  direct  strains  (longitudinal)  at  the  first  maximum  load  in 
compression.  As  for  elastic  analyses,  the  maximum  strain  occurs  at  the  hole  edge,  nearest 
to  the  'flanges',  at  half-thickness  of  the  plate.  The  finite  element  stress  history  (Von  Mises 
and  longitudinal  stresses)  at  the  position  of  maximum  stiain  is  shown  in  Figure  15  for  the 
CPLT  loading  cycle.  Figure  16  shows  the  variation  of  longitudinal  stress  with 
longitudinal  strain  at  the  maximum  strain  position  for  the  CPLT  loading  cycle.  Figure  17 
shows  the  residual  (longitudinal)  stresses  in  the  region  of  the  hole  at  the  completion  of 
the  first  CPLT  cycle  (end  of  Load  Step  (v)).  The  maximum  occurs  at  the  same  place  as  the 
maximum  total  strain,  as  shown  in  Figure  14.  Residual  stresses  at  the  end  of  the  second 
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CPLT  cycle  (end  of  Load  Step  (ix))  are  shown  in  Figure  18.  Comparison  of  maximum  and 
minimum  residual  stresses  after  the  first  and  second  CPLT  cycles  shows  an  increase  in  the 
tensile  stress  of  about  3%,  and  an  increase  in  compressive  stress  of  about  10%. 
Examination  of  Figures  17  and  18  shows  that  the  regions  of  tensile  and  compressive 
residual  stress  have  grown  from  the  first  to  the  second  cycle.  Figure  19  shows  the  residual 
Von  Mises  stress  distribution  along  the  half-width  of  the  specimen  at  the  end  of  the 
second  CPLT  cycle.  The  distribution  is  plotted  at  the  free  surface  of  the  specimen,  at  the 
mid-surface  and  at  Vi- thickness. 


5.  Conclusions 

A  Ramaswamy-Bodner-Stouffer  version  of  the  unified  constitutive  model  has  been 
successfully  implemented  into  the  ABAQUS  finite  element  package.  All  relevant 
combinations  of  numbers  of  direct  stresses  and  shear  stresses  have  been  tested, 
demonstrating  that  all  ABAQUS  (V5.6)  element  types  can  use  the  unified  constitutive 
model.  Results  from  a  suite  of  test  cases  agreed  with  answers  from  TREACLE,  the 
standalone  program  for  running  the  unified  constitutive  model.  TREACLE,  which 
outputs  stress  given  total  strain,  was  used  in  prior  work  to  validate  results  from  the 
PAFEC  implementation  of  the  model  [3],  verifying  that  the  ABAQUS  version  of  the 
model  is  consistent  with  the  PAFEC  version.  Therefore,  results  from  ABAQUS  and 
PAFEC  analyses  of  the  same  problem  can  be  compared  in  the  knowledge  that  material 
behaviour  is  modelled  identically.  Documentation  for  TREACLE  has  been  included  in 
this  report. 

Elastic  analyses  of  the  proposed  test  specimen  have  shown  that  the  specimen  should  be 
modelled  to  its  full  length,  with  the  load  applied  as  an  end  pressure,  and  that  the  contact 
pressure  exerted  by  the  test  machine  grips  can  be  ignored.  An  elastic  buckling  analysis 
predicted  that  the  correctly  aligned  and  fully  restrained  specimen  could  support  the 
proposed  test  load  of  about  1000  kN  without  excessive  lateral  displacements. 

The  ABAQUS  implementation  of  the  unified  constitutive  model  was  successfully  applied 
in  the  analysis  of  the  proposed  test  specimen.  This  showed  that  the  constitutive  model 
could  be  used  in  ABAQUS  for  real  problems  with  a  large  number  of  degrees-of-freedom. 
The  residual  stresses  predicted  under  simulated  CPLT  loading  verify  that  the  specimen  is 
suitable  for  experimental  testing.  It  is  recommended  that  an  experimental  program  be 
undertaken  on  the  specimen,  with  CPLT  loading  based  on  the  finite  element  analysis 
loading  scheme  (Section  4.2.3). 
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Appendix  A:  Unified  Constitutive  Model  Standalone 

Program  (TREACLE) 

A.1  TREACLE  Program 

TREACLE  is  a  standalone  program  which  'runs'  the  Ramaswamy-Bodner-Stouffer 
unified  constitutive  model,  in  the  form  described  in  this  report.  As  input  it  uses  material 
parameters  describing  the  constitutive  properties  of  the  material,  and  a  loading  history, 
which  is  typically  a  total  strain  rate  acting  over  a  specified  time.  The  input  parameters  are 
used  in  the  unified  constitutive  model  to  calculate  total  stress,  state  variables  and  other 
results  The  solution  process  is  controlled  by  a  series  of  commands  read  from  a  command 
file  the  format  of  which  is  based  on  an  ABAQUS  input  deck,  viz.,  command  lines 
beginning  with  *,  comment  lines  beginning  with  **,  and  data  lines.  User-selected 
variables  are  written  to  a  text  file  at  discrete  times  in  the  loading  history.  Finally,  plotting 
software  can  be  run  to  produce  a  graphical  representation  of  results  from  the  model. 

The  primary  purpose  of  TREACLE  is  to  aid  in  estimating  material  constants  for  use  in  a 
finite  element  method  implementation  of  the  unified  constitutive  model,  by  predicting 
values  of  stress/  strain  variables  which  can  readily  be  compared  to  experimental  results 
(eg  191).  The  PERTURB  function  in  TREACLE  allows  a  selected  parameter  to  take  a  range 
of  values  (with  other  parameters  held  constant),  to  assess  the  effect  of  varying  that 
parameter.  TREACLE  is  also  a  tool  for  validating  results  from  finite  element  method 
implementations  of  the  unified  constitutive  model.  For  instance,  the  Uniaxial  Stress-Strain 
Curve  in  Figure  1,  and  validation  solutions  in  Figures  3  and  5  were  obtained  by  running 
TREACLE  with  the  same  material  constants  and  total  strain  rate  as  the  finite  element  test 
cases.  Figure  A.l  contains  a  list  of  the  names  and  functions  of  the  routines  that  make  up 
TREACLE,  and  Figure  A.2  shows  the  calling  hierarchy  of  the  program. 
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Main  Program 

TREACLE  Displays  banner  and  opens  input  command  file 


Subroutines 

INITIAL 

INTERP 

PCHANGE 

DISPATCH 

RMATL 

SMATL 

INTTEMP 

RUNOUT 

OUTPUTI 

OUTPUTC 

OUTPUTR 

OUTPUTA 

PERTURB 

STORE 

RESTORE 

RUNSOLVE 

LOADUSE 

SOLVE 

LSODE 

DYDT 

CONLOAD 

EFFRAT 

EFFSTR 

CALCK2 

CALCWRK 

INTDY 

OUTPUTW 


Sets  initial  stresses  and  strains  to  zero 
Reads  command  lines 

Changes  values  of  material  constants,  as  required 

Interprets  command  lines 

Reads  material  properties  file 

Sets  variables  for  material  properties 

Interpolates  material  properties  with  respect  to  temperature 

Interprets  OUTPUT  command 

Initialises  output  files 

Closes  an  output  file 

Re-zones  an  output  file 

Appends  a  file  to  an  output  file 

Processes  PERTURB  command  by  looping  through  material  properties 

Stores  values  of  state  variables 

Restores  stored  values  of  state  variables 

Supplies  initial  conditions  to  the  solver 

Calculates  current  'load' 

Sets  parameters  for  the  Livermore  Solver 

Livermore  Solver  for  systems  of  Ordinary  Differential  Equations 

Calculates  state  variable  derivatives  for  LSODE 

Calculates  stress,  total  strain  and  elastic  strain 

Calculates  equivalent  inelastic  strain  rate 

Calculates  effective  stress 

Calculates  K2  for  inelastic  strain  rate  equations 

Calculates  inelastic  work  rate 

Interpolates  derivatives  of  state  variables  (part  of  Livermore  Solver) 
Writes  to  the  output  file 


Figure  A.l  FORTRAN  routines  in  unified  constitutive  model  standalone  program  TREACLE. 
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TREACLE 
|  INITIAL 
|  INTERP 

|  PCHANGE 

|  DISPATCH 
|  RMATL 
j  SMATL 
|  |  INTTEMP 

| RUNOUT 
j  |  OUTPUTI 

j  j  OUTPUTC 

|  |  OUTPUTR 

|  |  OUTPUTA 

j  PERTURB 
|  |  INTERP 

|  RESTORE 
|  STORE 
j  RUNSOLVE 

|  LOADUSE 
|  INTTEMP 
|  SOLVE 

|  LSODE* ** 
i  |  DYDT 

|  |  |  LOADUSE 

|  |  |  INTTEMP 

!  |  |  CONLOAD 

j  |  |  EFFRAT 

|  |  j  EFFSTR 

|  |  |  CALCK2 

j  j  |  EFFRAT 

j  |  |  CALCWRK 

j  LOADUSE 
|  INTTEMP 
j  CONLOAD 
|  INTDY" 
j  EFFRAT 
|  EFFSTR 
|  CALCWRK 
|  OUTPUTW 

Figure  A.  2  Sequence  of  routines  for  unified  constitutive  model  standalone  program  TREACLE 


*  LSODE  is  the  top-level  routine  of  the  Livermore  Solver  for  ODEs.  Its  call  structure  is  not  shown 
here,  except  for  DYDT,  which  is  a  user-supplied  subroutine. 

**  INTDY  is  part  of  the  Livermore  Solver,  but  here  it  is  called  independently 
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Figure  10.  Stress  mismatch  between  element  -predictions  at  common  nodes  (MPa)  -  initial 
coarse  finite  element  mesh  (1/8  model). 
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Figure  12.  Stress  mismatch  between  element  predictions  at  common  nodes  (MPa)  -  refined 
finite  element  mesh  (1/8  model). 
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Von  Mises  Stress  - Longitudinal  Stress  (Y-direction) 


Figure  15.  A  finite  element  stress  history  predicted  by  unified  constitutive  model  at 
critical  location  for  CPLT  loading. 


Longitudinal  Stress  vs  Longitudinal  Strain 


Figure  16.  Longitudinal  stress  vs  strain  predicted  by  unified  constitutive  model  at 
critical  location  for  CPLT  loading. 
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Figure  18.  Residual  longitudinal  stress  distribution  (<jyyr  MPa)  at  completion  of  2nd  CPLT 
cycle 
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(b)  Detail  near  local  minima  and  maxima 

Figure  19.  Residual  Von  Mises  stress  distribution  (MPa)  along  half-width  of  the  specimen  at 
completion  of  2nd  CPLT  cycle 
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